Initialisation et génération des données

## [1] "number of cores available = 1"
#Phi[1] ; eta = valeur de fin
#Phi[2] = valeur du noeud
#Phi[3] = echelle
m <- function(t, eta, phi) (phi[,1] + eta)/(1+exp((phi[,2]-t)/phi[,3]))
#=======================================#
param <- list(sigma2 = 0.05,
              rho2 = 0.1,
              mu = c(5,90,5),
              omega2 = c(0.5,0.1,0.01),
              #Survival data,
              nu2 = 0.5,
              a = 90,
              b = 50,
              alpha = 7,
              beta = 10)

F. <- length(param$omega2) #dimension de phi
#=======================================#
t <- seq(60,120, length.out = 10) #value of times

# --- longitudinal data --- #
G = 40 ; ng = 8 #nombre de groupe et d'individu par groupe
n <- G*ng*length(t) ; N <- G*ng  #nombre total de data, et #nombre total d'individu

dt_NLME <- NLME_obs(G, ng, t, param, m)
Y <- dt_NLME$obs
dt_SF <- SF_obs(dt_NLME, param, m)

Statistique exhaustive

\[\log f(Y, Z ; \theta) = \langle \Phi(\theta) ; S(Y,Z) \rangle - \psi(\theta)\]

sigma2_a <- sigma2_b <- sigma2_alpha <- 0.1
#Petite fonction pour retourner rapidement l'appel Phi[attr(Phi, i)] où i est '1', '2', ...,
`%a%` <- function(x,var){
  if(length(var)== 1) return(x[ attr(x,as.character(var)) ])
  lapply(var, function(v) x%a%v) %>% unlist
}


Phi <- fct_vector(function(sigma2, rho2, mu, omega2, nu2, a, b, alpha, beta) {
  c(- n/(2*sigma2), -N/(2*rho2),               #1, 2
        - G/(2*omega2), G*mu/omega2,               #3, 4
        - G/(2*nu2),                               #5
          a/sigma2_a,         -1/(2*sigma2_a),     #6, 7
          b/sigma2_b        , -1/(2*sigma2_b),     #8, 9
          alpha/sigma2_alpha, -1/(2*sigma2_alpha), #10, 11
          beta ) },#12
  dim = c(1,1,F., F., rep(1,8)) )

zeta <- function(beta, eta, phi, gamma, a,b, alpha)
{
  lbd <- function(t,g) t^{b-1} * exp(alpha*m(t, eta[g], phi[g,]))
  P1 <- 1:nrow(dt_SF) %>% sapply(function(i) integrate(function(t) lbd(t,dt_SF$gen[i]) , 0, dt_SF$obs[i])$value )
  
  P2 <- 1:nrow(dt_SF) %>% sapply(function(i) exp(beta*dt_SF$U[i] + gamma[dt_SF$gen[i]] )  )
  
  beta*sum(dt_SF$U) -  b*a^-b * sum(P2*P1)
}
zeta.der <- function(beta, eta, phi, gamma, a, b, alpha)
{
  lbd <- function(t,g) t^{b-1} * exp(alpha*m(t, eta[g], matrix(phi[g,], nrow = 1)))
  P1 <- 1:nrow(dt_SF) %>% sapply(function(i) integrate(function(t) lbd(t,dt_SF$gen[i]) , 0, dt_SF$obs[i])$value )
  
  P2 <- 1:nrow(dt_SF) %>% sapply(function(i) dt_SF$U[i]*exp(beta*dt_SF$U[i] + gamma[dt_SF$gen[i]] )  )
  
  sum(dt_SF$U) -  b*a^-b * sum(P2*P1)
}

S <- fct_vector(function(eta, phi, gamma, a, b, alpha) mean((Y - get_obs(m, dt_NLME, eta = eta, phi = phi) )^2 ),
                function(eta, phi, gamma, a, b, alpha) mean(eta^2),           #2
                function(eta, phi, gamma, a, b, alpha) apply(phi^2, 2, mean), #3
                function(eta, phi, gamma, a, b, alpha) apply(phi, 2, mean),   #4
                function(eta, phi, gamma, a, b, alpha) mean(gamma^2),         #5
                function(eta, phi, gamma, a, b, alpha) a,                     #6
                function(eta, phi, gamma, a, b, alpha) a^2,                   #7
                function(eta, phi, gamma, a, b, alpha) b,                     #8
                function(eta, phi, gamma, a, b, alpha) b^2,                   #9
                function(eta, phi, gamma, a, b, alpha) alpha,                 #10
                function(eta, phi, gamma, a, b, alpha) alpha^2,               #11
                function(eta, phi, gamma, a, b, alpha)                        #12
                  uniroot(function(beta)zeta.der(beta, eta, phi, gamma, a,b,alpha), lower = -100, upper = 100)$root,
                dim = c(1,1,F., F., rep(1,8)) )

Metropolis Hastings

Vraisemblance

loglik.phi <- function(x, eta, Phi){
  id <- c(1,3,4)
  sum( Phi%a%id * S(eta = eta, phi = x, i = id) )
}
loglik.eta <- function(x, phi, Phi){ 
  id <- c(1,2)
  sum( Phi%a%id * S(eta = x, phi = phi, i = id) )
}

loglik.gamma <- function(x, eta, phi, a, b, alpha, Phi){
  Phi%a%5  * sum(x^2) + ng*sum(x) +
    zeta(Phi%a%12, eta, phi, x, a, b, alpha) #Correction fait à la va vite
}

loglik.a <- function(x, b, Phi) sum( Phi%a%6:7 * S(a = x, b = b, i = 6:7) ) + G*ng*log(b*x^-b)
loglik.b <- function(x, Phi)    sum( Phi%a%8:9 * S(b = x, i = 8:9) ) + (x-1)*sum(log(t))
loglik.alpha <- function(x, Phi)sum( Phi%a%10:11 * S(alpha = x, i = 10:11) )

SAEM

initialisation

M <- 1 #nombre de simulation

# ---  Initialisation des paramètres --- #
para <- param %>% sapply(function(x) x + runif(1))
para$rho2 = 0.2 ; para$mu <- c(6,3,1) ; para$omega2 <- rep(.1,3)

# --- Initialisation des chaines MC : Z_0 --- #list(attr(dt_NLME,'eta')),#list(attr(dt_NLME,'phi')),#
Z <- list(eta = 1:M %>% lapply(function(i) rnorm(G*ng, 0, para$rho2)  %>% matrix(ncol = 1) ),
          phi = 1:M %>% lapply(function(i) matrix(rnorm(F.*G, para$mu, para$omega2), nrow = F.) %>% t ),
          
          gamma = list(attr(dt_SF,'gamma')),#1:M %>% lapply(function(i) matrix(rnorm(G, 0, para$nu2), ncol = 1) ),
          
          a = list(matrix(para$a)),
          b = list(matrix(para$b)),
          alpha = list(matrix(para$alpha)) )

Etape simulation et maximisation du SEAM

sim <- function(niter, Phih, eta, phi, gamma, a, b, alpha)
{
  M <- length(phi)

  eta <- 1:M %>% lapply( function(i)
    MH_High_Dim_para_future(niter, eta[[i]], sd = .036,                 loglik.eta, phi[[i]], Phih, cores = ncores ))

  phi <- 1:M %>% lapply( function(i)
    MH_High_Dim_para_future(niter, phi[[i]], sd = c(.028, .036, .021), loglik.phi, eta[[i]], Phih, cores = ncores ))
  
  # gamma <- 1:M %>% lapply( function(i)
  #   MH_High_Dim_para_future(niter, gamma[[i]], sd = .03, loglik.gamma, 
  #                    eta[[i]], phi[[i]], a[[i]], b[[i]], alpha[[i]], 
  #                    Phih, cores = 1 ))

  a <- 1:M %>% lapply( function(i) MH_High_Dim_para_future(niter, a[[i]],     sd = .02, loglik.a,b[[i]],  Phih, cores = 1 ))
  b <- 1:M %>% lapply( function(i) MH_High_Dim_para_future(niter, b[[i]],     sd = .02, loglik.b,         Phih, cores = 1 ))
  alpha <- 1:M %>% lapply( function(i) MH_High_Dim_para_future(niter, alpha[[i]], sd = .02, loglik.alpha, Phih, cores = 1 ))

  list(eta = eta, phi = phi, gamma = gamma, a = a, b = b, alpha = alpha)
}

maxi <- function(S)
{
  list(sigma2 = S%a%1,
       rho2 =   S%a%2,
       mu =     S%a%4,
       omega2 = S%a%3 - (S%a%4)^2,
       nu2 =    S%a%5,
       a =      S%a%6,  b =    S%a%8,
       alpha =  S%a%11, beta = S%a%12 )
}

Resultats

sigma2 rho2 mu1 mu2 mu3 omega21 omega22 omega23 nu2 a b alpha beta
Oracle 0.051138 0.0995415 4.863004 90.04573 5.032278 0.4684432 0.0986736 0.0075891 0.4557114 90.00000 50.00000 49.000000 9.96179
Initialisation 0.414923 0.2000000 6.000000 3.00000 1.000000 0.1000000 0.1000000 0.1000000 1.1002253 90.94465 50.68087 7.226237 10.73202 
niter <- 100
MH.iter <- 10
res <- SAEM(niter, MH.iter, para, Phi, S, Z, sim, maxi, eps = 1e-3, verbatim = 2)
saveRDS(res, 'saem.rds')
gg <- plot(res)
## [1] "SAEM execution time = 06min 56sec"
Résultat de l’algo SAEM-MCMC
sigma2 rho2 mu.1 mu.2 mu.3 omega2.1 omega2.2 omega2.3 nu2 a b alpha beta
Valeur réelle 0.0500 0.1000 5.0000 90.0000 5.000 0.5000 0.1000 0.0100 0.5000 90.0000 50.0000 7.0000 10.0000
Valeur estimée 0.1258 4.6226 3.0392 3.0394 0.965 0.1165 0.0031 0.0237 0.4557 83.7281 55.1365 231.8184 -2.0753
Rrmse 1.5162 45.2257 0.3922 0.9662 0.807 0.7669 0.9690 1.3671 0.0886 0.0697 0.1027 32.1169 1.2075